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ABSTRACT 


Computational fluid dynamics is being used increasingly to predict flows for 
aerospace propulsion applications, yet there is still a need for an easy to use, computa- 
tionally inexpensive turbulence model capable of accurately predicting a wide range 
of turbulent flows. The Baldwin-Lomax model is the most widely used algebraic 
model, even though it has known difficulties calculating flows with strong adverse 
pressure gradients and large regions of separation. The modified mixing length model 
(MML) was developed specifically to handle the separation which occurs on airfoils 
and has given significantly better results than the Baldwin-Lomax model. The success 
of these calculations warrants further evaluation and development of MML. 

The objective of this work was to evaluate the performance of MML for zero 
and adverse pressure gradient flows, and modify it as needed. The Proteus Navier- 
Stokes code was used for this study and all results were compared with experimental 
data and with calculations made using the Baldwin-Lomax algebraic model, which is 
currently available in Proteus. 

The MML model was first evaluated for zero pressure gradient flow over a flat 
plate, then modified to produce the proper boundary layer growth. Additional modifi- 
cations, based on experimental data for three adverse pressure gradient flows, were 
also implemented. The adapted model, called MMLPG (modified mixing length 
model for pressure gradient flows), was then evaluated for a typical propulsion flow 
problem, flow through a transonic diffuser. Three cases were examined: flow with no 
shock, a weak shock and a strong shock. 



The results of these calculations indicate that the objectives of this study have 
been met. Overall, MMLPG is capable of accurately predicting the adverse pressure 
gradient flows examined in this study, giving generally better agreement with experi- 
mental data than the Baldwin-Lomax model. 
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CHAPTER I 


INTRODUCTION 

1.1 Motivation and Objectives 

Computational Fluid Dynamics (CFD) is a valuable tool for calculating the 
turbulent flow fields that occur in engineering fluid flow problems. Some of the 
characteristics of turbulent flow include random fluctuations in fluid properties, the 
enhancement of mixing, diffusion and dissipation, and the presence of eddies of 
various sizes. Turbulent flow is, therefore, very difficult to predict theoretically. 
Experiments provide much useful information about turbulent flow fields but are 
costly and time consuming, so CFD is being used increasingly to reduce or optimize 
the amount of experimental testing which must be done. 

Most CFD codes solve the equations of conservation of mass, momentum 
(Navier-Stokes) and energy and, in principle, completely describe the details of turbu- 
lent flow. However, except for very simple problems, these equations cannot be 
solved exactly due to the limited capabilities of computational resources. Most 
engineering problems are primarily concerned with mean fluid properties and not with 
the details of the turbulent fluctuations; the mean properties can' therefore be 
computed using the Reynolds-averaged form of the Navier-Stokes equations. 1 In 
Reynolds averaging, the conservation equations are averaged over a time scale that is 
large compared to the largest time scale of the fluctuating motion. 1,2,3 The averaging 
procedure introduces new terms which represent the turbulent transport of mean 
momentum, heat and mass. The resulting averaged equations are not closed and 
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empirical information, in the form of a turbulence model, must be used to close the 
system. 

A turbulence model is a mathematical model consisting of an equation or set 
of equations which determines the turbulent transport terms in the mean flow 
equations and hence closes the system of equations. 1 Turbulence models give an 
approximate description of the flow by describing the overall effect of turbulence on 
the mean flow, rather than describing the details of the turbulent motion. Since turbu- 
lent transport processes depend on factors such as geometry, swirl effects and 
buoyancy, turbulence models, which are usually developed based on hypotheses about 
a certain flow or range of flows, usually have a limited range of applicability. 
Typically, a model which is complex and consists of a large number of equations is 
difficult to use and is computationally expensive. Often this increase in “cost” is not 
proportional to the improvements in the computation. 

For most engineering applications, a turbulence model should be easy to 
implement, computationally inexpensive and applicable to a wide range of flows. 
Algebraic turbulence models, also called zero-equation models, are simple and 
inexpensive, however they generally have only a narrow range of applicability. The 
most widely used algebraic model, the Baldwin-Lomax model (BLM), 4 fits this 
description, but it is known to have difficulties calculating adverse pressure gradient 
and separated flows, 5 " 12 the regime it was designed to handle. 

In 1989, the modified mixing length model (MML) was developed and used to 
calculate separated flows over airfoils, flowfields that BLM was unable to accurately 
predict. 5 It is based on Prandtl’s mixing length hypothesis 3 and uses a mixing length 
that is dependent on the local wall shear stress. The objective of this work is to 
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continue the development of MML to expand its range of applicability to include 
boundary layer flows with adverse pressure gradients. 

1.2 Overview 

Chapter II gives some background information on the Proteus Navier-Stokes 
code/ which was used to make all of the calculations in this work. It also describes the 
implementation of turbulence into the governing equations and describes problems 
encountered with BLM, the current algebraic turbulence model in Proteus. Chapter II 
also describes the original formulation of MML. Chapter III reports calculations 
made with MML for zero pressure gradient flow over a flat plate, and then describes 
the modifications made to improve these results for both zero and adverse pressure 
gradient flows. The resulting modified version of MML is called MMLPG. Chapter 
IV compares MMLPG and BLM for three transonic diffuser flow test cases: flow with 
a weak shock, strong shock, and no shock. Chapter V contains a summary of this 
work and a discussion of the conclusions drawn. 



CHAPTER II 


BACKGROUND 

2.1 The Proteus Navier-Stokes Code 

The Proteus Navier-Stokes code, 13, 14 developed at the NASA Lewis 
Research Center, is a user-oriented, foil Navier-Stokes code for aerospace propulsion 
applications. Proteus solves the Reynolds-averaged, unsteady, compressible Navier- 
Stokes equations in strong conservation law form. Two separate versions of the code 
exist: one for two-dimensional plane or axisymmetric flow, and one for three-dimen- 
sional flow. A primary objective of the Proteus effort was to make the code easy to 
use and modify. Therefore, code readability, modularity and documentation were 
emphasized, rendering the code ideal for the insertion and development of a new 
turbulence model. 

The governing equations in Proteus are written in Cartesian coordinates and 
then transformed to a nonorthogonal, body-fitted system (see Appendix l). 13 They 
are solved by marching in time using a folly-coupled alternating direction implicit 
solution procedure with generalized first or second order time differencing. 15, 16 The 
boundary conditions are also treated implicitly and can be steady Or unsteady. All 
terms, including diffusion terms, are linearized to second order using Taylor series 
expansions. The two turbulence models originally available in Proteus are the 
Baldwin-Lomax algebraic model 4 and the Chien k-e two-equation model. 17 

In addition to solving the full, time-averaged Navier-Stokes equations, 
Proteus includes options to solve the thin-layer and Euler equations, and to eliminate 
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the energy equation by assuming constant stagnation enthalpy. Artificial viscosity is 
used to minimize the odd-even decoupling resulting from the use of central spatial 
differencing for the convective terms, and to control pre- and post-shock oscillations 
in supersonic flow. 13 Two artificial viscosity models are available: a combination 
implicit/explicit constant coefficient model, 18 and an explicit nonlinear coefficient 
model designed specifically for flows with shock waves. 19 The artificial viscosity is 
discussed in more detail in Appendix 2. At the NASA Lewis Research Center the 
code is typically run either on the CRAY X-MP or CRAY Y-MP computer, and is 
highly vectorized. For all calculations made herein, the two-dimensional/axisymmet- 
ric version of the code was run on the CRAY Y-MP computer. 

2.2 Algebraic Turbulence Modeling and the Baldwin-Lomax Model 

Accurate modeling of turbulence is essential to the computation of complex 
propulsion flow fields. Several types of turbulence models are available, ranging 
from zero-equation algebraic models to multi-equation Reynolds-stress models. 
Algebraic models are the most algorithmically simple and computationally inexpen- 
sive models and were therefore chosen as the focus of this effort. 

Proteus, along with the majority of Navier-Stokes codes, uses the Boussinesq 
assumption, 3 which states that the turbulent stresses behave like the molecular viscous 
stresses and therefore are proportional to the mean velocity gradient. The resulting 
total shear stress for a two-dimensional flow is given by 13 

* = < 21 > 

The effective viscosity is defined as jt eff = jx + |i t , where |X is the molecular viscosity 
and p t is the turbulent, or “eddy” viscosity. The same analogy applies to the heat flux 
and the normal stresses, which are both defined in Appendix 1, such that an effective 
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second coefficient of viscosity is defined as A, eff = X + X t and an effective thermal 
conductivity coefficient is defined as k eff = k + k t . 

Most algebraic turbulence models are based on Prandtl’s mixing length 
hypothesis which builds on the Boussinesq assumption. 1 Prandtl made an analogy 
between molecular motion and turbulent flow. In molecular motion, the molecular 
viscosity is proportional to the average velocity and the mean free path of the 
molecules. In turbulent flow, Prandtl assumed that the turbulent viscosity is propor- 
tional to the characteristic velocity of the fluctuating motion and to a typical length, 
called the “mixing length”, of this motion. In other words, 

(x t = pv t / (2.2) 

where v, is the turbulent velocity scale and the mixing length, l, is the transverse 
distance over which fluid particles maintain their original momentum. Prandtl further 
assumed that the turbulent velocity scale is equal to the mixing length times the veloc- 
ity gradient so that 





(2.3) 


The quantity l 
primary flow d 
tion. 


5u 

dy 


is the velocity scale, where u is the component of velocity in the 


irection and y is the coordinate perpendicular to the primary flow direc- 


The current algebraic turbulence model in Proteus, the Baldwin-Lomax model 
(BLM), is given in Appendix 3. It is the most well-known and widely used algebraic 
turbulence model. An extension of the Cebeci-Smith model, 20 which requires knowl- 
edge of the outer edge of the boundary layer, the Baldwin-Lomax model was devel- 
oped to handle separated flows while avoiding the necessity of finding this outer edge. 
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Several references report problems with BLM in regions of strong pressure 
gradient and in flows with large regions of separation. Yu 6 reports problems calculat- 
ing surface pressures on the outboard wing region of a wing-body configuration when 
the angle of attack is high and separation occurs. Visbal and Knight' report that the 
BLM outer formulation is unsuitable for separated supersonic flow and also is unable 
to predict the recovery of the boundary layer downstream of reattachment. Degani & 
Schiff 8 report that BLM is unsuitable in regions of cross flow separation due to 
ambiguities in computing the outer length scale. Menter 9 reports that BLM underesti- 
mates the displacement thickness for the increasingly adverse pressure gradient flow 
of Samuel & Joubert 21 and that it also gives an incorrect prediction of the adverse 
pressure gradient flow of Driver. 22 Potapczuk 5 reports problems with BLM in 
predicting the separation and unsteady behavior on airfoils with and without leading 
edge ice accretion. Stock and Haase 10 report that BLM does not predict the correct 
trends for p t or the Reynolds stresses in adverse pressure gradient and separated 
flows. 

There are several reasons why BLM has problems computing flows with large 
pressure gradients and large regions of separation. The primary difficulty occurs in 
finding the maximum of the F(y) function (defined Appendix 3), which has two or 
more peaks in regions of separated flow. This behavior is shown in figure 1 taken 
from reference 5. As the relative magnitudes of the local maxima change, y max 
(defined in Appendix 3), may suddenly jump, producing unrealistic discontinuities in 
the turbulent viscosity. Selection of the global maxima often results in a gross over- 
prediction of the turbulent viscosity. Some authors 7, 23 have found that choosing the 
outermost peak produces better results, while others have elected to use the innermost 
peak. 8 To account for upstream turbulence history effects, Visbal and Knight 7 used 
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BLM with relaxation. They also found that the BLM constants C cp and C Kleb (see 
Appendix 3) should vary with Mach number. Sakowski et al 11 expounded upon this 
by finding a relation for C cp as a function of Mach number and pressure gradient, but 
encountered problems caused by the vanishing of the Van Driest factor when x w , the 
local shear stress at the wall, approaches zero. To remedy this, they used the local 
shear in place of x w . Launder and Pridden 24 report several modifications to the Van 
Driest factor 2 for pressure gradient flows, many of which incorporate the local shear 
stress instead of x w . 

A simpler, yet effective, approach was used by Potapczuk, 5 who developed a 
modified mixing length (MML) model that does not require a boundary layer thick- 
ness, but also avoids all problems associated with the determination of a maximum 
F(y). This model is described in detail in the following section. 

2.3 The Modified Mixing Length Turbulence Model 

The modified mixing length (MML) model was developed by Potapczuk 5 to 
fill the need for an algebraic model to handle turbulent flow with large separated 
regions. The particular problem of interest was an airfoil at angle of attack with and 
without leading edge ice accretions. Previous calculations made with BLM gave poor 
results and the source of the problem was the function F(y), which had multiple peaks 
for this flow case. 

The MML model avoids the need to seek a maximum of some ad hoc function. 
In accordance with Prandtl’s mixing length theory, the MML model determines the 
mixing length using the wall shear stress and the normal distance from the wall, with 
the maximum mixing length capped off at a given value. Thus, it is a two layer 
model, such that the length scale depends on conditions near the surface and remains 
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constant in the separated region. This assumption is valid since there is no substantial 
enhancement of turbulence in separated regions. The turbulent viscosity is given by 

p t = pflcoj (2.4) 

where the velocity gradient in equation (2.3) has been replaced by the vorticity magni- 
tude, }co] . Figure 2, taken from reference 20, shows the behavior of the mixing length 
in a turbulent boundary layer. Several empirical formulas are available to evaluate the 
inner region, 20, 24 which consists of the viscous sublayer and the overlap layer. The 
MML model uses the van Driest formulation, 20 which is given by 


= Ky ^1 - < 


-A 


(2.5) 


where A + = 26 and the value of k, the von Karman constant, is 0.4. The quantity y + 
is defined as 



where y is the shear length scale. 


( 2 . 6 ) 


JpK\ 


(2.7) 


For y + £ 5 A + (but still in the “inner” region), the mixing length is approximated by 
Ky; this is the original Prandtl theory, and is consistent with the well-known logarith- 
mic profile. In the outer region of the boundary layer, the outer mixing length behaves 
according to 


Pouter = constant xy* 


( 2 . 8 ) 
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The MML model uses a blending function to give a smooth transition between 
the inner and outer layers and is given by 

/(y) = K^y* (l - (1 - £-) ] (l - e^j, y* < C, (2.9) 

i(y) = K^iy*, y + >C, (2.10) 

In this formulation, C,y is the distance above the surface at which l saturates, and C 2 
controls the curvature of the blending region. See figure 3 for a typical MML model 
mixing length profile. 

Calculations made by Potapczuk with the MML model showed improvements 
over the BLM calculations for the prediction of the separated region, the maximum 
lift coefficient and vortex shedding frequencies. Since the MML model was devel- 
oped to solve the specific problem of flow over airfoils, a comprehensive evaluation 
of the model for more general flowfields was not a part of that study. The objective of 
the present study is to evaluate the MML model for general zero pressure gradient and 
adverse pressure gradient turbulent boundary layer flows and examine possible 
modifications to improve the performance of the model. 




CHAPTER IH 


EVALUATION AND MODIFICATION OF MML 


The MML model, described in Chapter II, was modified so that it could better 
calculate turbulent boundary layer flows with zero and adverse pressure gradients. 
The first step in this process was to optimize the computation of the wall shear stress 
used in MML. Next, MML was used to calculate a turbulent boundary layer with zero 
pressure gradient. These results exhibited poor agreement with experimental data, so 
modifications were made to MML to remedy this problem. Then modifications corre- 
sponding to two adverse pressure gradient flows of Bradshaw 25 were successfully 
incorporated into MML. Finally, all of the modifications were combined into one 
general model, called MMLPG. 

3.1 Optimization of Shear Stress Estimate 

Since the MML turbulence model is a function of the wall shear stress, it is 
important to accurately calculate this quantity. The wall shear stress is given by 

X. = OS> 


The molecular viscosity, p, is a very small quantity compared to the velocity gradient, 

(i^-) , which is strongly dependent on several factors such as the finite difference 

oy w 

scheme used, the grid spacing and the numerical features of the code. It is important 

to minimize the sensitivity of these factors because small changes in the estimate of 
3u 

(• 5 —) may actually produce large changes in T w . A more global approach is to use a 
oy w 
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parabolic extrapolation of x w , using the shear stress at two interior grid points, and the 
momentum equation in the streamwise direction, which reduces to 


W„ - 


at the wall. Using a parabola to define x. 


x = ay +by + c 


where. 


(3.2) 


(3.3) 


a 


x I -x 2 -b(y 1 -y 2 ) 

y?-yi 


(3.4) 


b = 



(3.5) 


c = i[x, + x 2 -a(y? + y|) -b(y, + y 2 )] (3.6) 

Here, the subscripts 1 and 2 denote interior points as depicted in figure 4. Note that 
x w =c, since y=0 at the wall. Also note that the shear stress at interior points is defined 
by 


,du , 3v. 
X “ ^ {f( dy dx ) 


(3.7) 


The parabolic extrapolation in equations (3.3) through (3.6) gives a reliable value for 

3u 

x w and avoids problems that could arise from sensitivity of the (^-) estimate. In 

W 5u 3u ^ w 

fact, (^-) can be found from (^— ) = x w /|i. 

oy w w 
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Modifications were made to the Proteus code to calculate the shear stress 
profile in the boundary layer. These modifications made use of the generalized grid 
transformations in Proteus such that 


9u , 3u t 3v . _0v. 
T - ^eff + + 


(3.8) 


where £ y = ri x = 0 for an orthogonal grid. 13 The derivatives in the above equation 
are calculated in Proteus using 3-point, second-order central differencing. To see if 
higher order differencing would improve the calculation, 5-point, fourth-order central 
differencing 26 was also used to calculate the velocity gradients. 

The test case of incompressible, zero pressure gradient, turbulent flow over a 
flat plate, as shown in figure 5, was used to evaluate the shear stress calculations. The 
grid, shown in figure 6, had 51 points in both the streamwise and normal directions 
and had grid points clustered at the wall to resolve the boundary layer and at the 
upstream boundary to resolve the imposed boundary condition. In addition, it was 
evaluated to insure grid indepence for zero pressure gradient flow. The reference 
velocity, temperature, pressure and length used in Proteus were 33.53 m/s, 288.3 K, 
101.3 kPa and 1.98 m, respectively. At the upstream boundary, the velocity profile, 
which was computed using the correlation of Musker, 27 was held fixed. The flow was 
computed using both MML and BLM, using both higher and lower order differencing 
of the velocity gradients in the shear stress computation. The MML constants were 
chosen as Ci=3000 and C2=5, which were found to give good results at Re x =7xl0 6 . 
Both turbulence models produced good agreement with experimental velocity-defect 
profiles, 28 as shown in figure 7, which shows calculations at Re x =7xl0 6 made with 
the lower order differencing of the velocity gradients. The quantity u T in figure 7 is 


the shear velocity, given by u t = J(\x w \/p) . The accuracy of the finite differencing 
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Figure 7. Velocity-defect profiles for zero pressure gradient flat plate flow, 
Re x =7xl0 6 


. 20 

produced no noticeable improvement in the velocity profiles, but there were slight 
differences in the shear stress profile very close to the wall, which can be observed in 
the plots of figure 8. The higher-order differencing produced a smoother shear stress 
profile near the wall, and thus subsequent calculations will make use of the higher- 
order calculation of the velocity gradients. 

Near separated regions of flow, x w approaches zero which will cause y and 
thus the mixing length to become infinitely large. To avoid this problem, the follow- 
ing local average used by Potapczuk 5 was incorporated: 

M = 0.1|VJ+0.2|V,| +0.4|tJ +0.2|t,J +0.1|t, J (3.9) 

The subscripts in equation (3.9) refer to grid points in the streamwise direction along 
the wall. 

3.2 Evaluation and Modification for Zero Pressure Gradient Flows 

In reference 5, a series of cases were run for flow over a NACA0012 airfoil at 
conditions near stall using both BLM and MML. The MML constants Ci=2000 and 
C2=5 were chosen based on correlations with experimental data. The Baldwin-Lomax 
model tended to suppress the trailing-edge separation, which occurs on the top surface 
of the airfoil, by over-predicting p. t throughout the separated region. On the other 
hand, MML predicted high values of ji t only near the separation point, thus allowing 
the reverse flow to develop downstream. In the current study, MML was evaluated for 
turbulent flow over a flat plate at zero pressure gradient, as described below. 

In the preliminary analysis presented in section 3.1, the constants Cj = 3000 
and C 2 = 5 were found to give good agreement for Re x = 7xl0 6 . At other locations on 
the plate, i.e., at other Reynolds numbers, the BLM velocity-defect profiles correctly 
exhibit similarity but the MML profiles do not, as shown in figure 9. In a turbulent 
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y/8 


(a) MML 



y/5 


(b) BLM 


Figure 9. Velocity-defect profiles for zero pressure gradient flat plate flow at three 
Reynolds numbers. 
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flow over a flat plate, the boundary layer thickness increases with increasing x- 
distance along the plate. To accurately model this flow, the turbulent length scale 
must also increase proportionately with the boundary layer thickness. In MML, the 
outer length scale, as given in equation (2.10), is equal to a constant times the shear 
length scale, y . The increase in y with x-distance is negligible, resulting in an 
essentially constant value of the outer length scale for all Reynolds numbers. Plots of 
p t , as shown in figure 10, illustrate that the turbulent viscosity profiles calculated with 
MML are nearly the same at all Reynold’s numbers, but the BLM |X t profiles increase 
with increasing Reynolds number. Though MML produced the correct length scales 
for an airfoil near stall, 5 modifications are needed to make it applicable to general 
boundary layer flows. 

In order to make MML applicable over a range of Reynolds numbers, the 
optimal saturation lengths, or Cj values, were found at several Reynolds numbers. 
The following simplified formulas were used to calculate the inner and outer mixing 
lengths: 

t = K y + 1^1 — e A+ J, 

/ + cp = kC, C, £y + (3.11) 

Here, t is the nondimensional form of the mixing length, equivalent to l/y*, and the 
outer length scale, / + ctp , is simply the inner length scale evaluated at y + = C, . From 
these results, / + cip was found as a function of the skin friction, Cf, giving 


C,<y + 


(3.10) 


/ + ctp = 1 860 - 6.20 xl0 5 c f 


(3.12) 






25 


The velocity-defect profiles of figure 11 show that equations (3.10) and (3.12) with 
/ + = min (/ + , ? ct p ) , allow the mixing length to grow proportionately with the bound- 
ary layer thickness. The modified MML is better than BLM at predicting the local 
skin friction coefficient, Cf, as shown in figure 12; the wiggles at the upstream bound- 
ary are a result of the imposed upstream boundary condition. 

3.3 Modifications for Adverse Pressure Gradient Flows 

Two equilibrium pressure gradient flows of Bradshaw 25 were used to modify 
MML for adverse pressure gradients effects. Equilibrium turbulent flows are flows 
which have a constant value of Clauser’s equilibrium parameter, 

0 = (3.13) 

x w 3x 

In addition, they correspond to a power-law velocity profile distribution, u e x\ 
where the magnitude of the- exponent, a, indicates the strength of the pressure gradi- 
ent They also exhibit similarity when plotted in velocity-defect coordinates. Three 
flows were examined in the experimental study of Bradshaw; 25 these were flows with 
zero, mild, and strong pressure gradients. The corresponding values of the exponent, 
a, are 0, -0.15, and -0.255, respectively; the corresponding values of 0 are 0, 1 and 5. 

The modifications to the turbulence model are based on the trends exhibited in 
the mixing length at the three pressure gradients as shown in figure 13 taken from 
Bradshaw 25 . Note that for all three pressure gradients, the maximum mixing length is 
approximately 0.088, the saturation distance from the wall is approximately 0.48, and 
the slope of the curves near the wall increases with the strength of the pressure gradi- 
ent. These three features were used to develop the following model: 
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Figure 11. Velocity-defect for zero pressure gradient flat plate flow calculated using 
the modified MML of section 3.2. 
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Re e 


Figure 12. Local skin friction coefficient for zero pressure gradient flat plate flow; 
MML of section 3.2 used. 
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(3.14) 


/ = kG. 


G, 


y + > G, 


(3.15) 


A new parameter, G 3 , has been introduced and the constants and C 2 in the original 
MML model have been replaced by the functions Gj and G 2 , where 


G, = 0.4G 4 

(3.16) 

G 2 = 5G 3 k 

(3.17) 


Here, G 4 is essentially a nondimensional boundary layer thickness which is a function 
of (3 and Cf, and G 3 controls the slope of mixing length curve and is a function of {3. 
The following correlation was assumed for G 4 : 

G 4 = G 5 + G 6 c f (3.18) 

Separate values of the parameters G 3 , G 5 , and Gg, corresponding to each of the three 
pressure gradients, were found and are given in table 1. This results in essentially 
three separate models, one for each pressure gradient, depending on which set of 
parameters is used. The results of these -modifications are compared with Baldwin- 
Lomax calculations in the velocity-defect plots of figures 14 through 16. (Note: The 
Baldwin-Lomax results for the zero pressure gradient case are given in figure 9.) The 
reference conditions used are the same as those given in section 3.1. For the two 
adverse pressure gradient cases, the turbulent velocity profiles at the upstream bound- 
ary were computed using a cubic spline fit of the Bradshaw experimental data and 
held fixed; the appropriate pressure gradient was imposed at the freestream boundary. 
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Table 1. Parameters used in pressure gradient modifications. 


Pressure 

Gradient 

Strength 

P 

g 3 

g 5 

g 6 

zero 

0 

1.00 

23,300 

-7.75xl0 6 

mild 

1 

1.25 

30,100 

-1.16xl0 7 

strong 

5 

1.53 

33,800 

-2.09x1 0 7 



Figure 14. Velocity-defect for zero pressure gradient flow calculated using the 
modified MML of section 3.3. 
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Both cases were computed using the grid of figure 6, but for the strong pressure gradi- 
ent case, the number of grid points in the vertical direction was increased to 101. 

3.4 Final Model 

The final step in developing this turbulence model was to combine all of the 
above modifications to get one general turbulence model. To accomplish this, the 
following correlations were developed for the parameters G3, G5 and Gg: 


G 3 = 1.0 P<0.0 (3.18a) 

G 3 = 1.0 + 0.307(5- 0.0391 p 2 0.0<P<5.34 (3.18b) 

G 3 = 1.52 P > 5.34 (3.18c) 

G s = 23, 300 P<0.0 (3.18d) 

G 5 = 23, 300 + 8560p - 1230P 2 0.0 < p < 5.34 (3.18e) 

G s = 33,900 P > 5.34 (3.18f) 

G 6 = -7.75xl0 6 p<0.0 (3.18g) 

G 6 = - 7.75xl0 6 - 4.51 x 10 6 + 386, OOOp 2 0.0 < p < 5.34 (3.18h) 

G 6 = -20, 900 p > 5.34 (3.18i) 


The available experimental data is limited to only the three values of P which are in 
the range 0< P< 5.34, and the quadratic correlations of equations (3.18b), (3.18e) and 
(3.18h) are based on this limited data. For P<0, the values in (3.18a), (3.18d) and 
(3.18g), were obtained by evaluating the quadratic equations at P=0. Similarly, for 
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£>5.34, the values in (3.18c), (3.180 and (3.18i) were obtained by evaluating the 
quadratic equations at {3=5.34. 

Since (3, which is defined in equation (3.13), is a function of the displacement 
thickness, 8i, a correlation was also developed to avoid the problem of calculating Sj 
directly and thus having to define the edge of the boundary layer: 

8, = (G 7 + G 8 c f ) y* (3.19) 

The parameters G7 and Gg were defined in a manner similar to G3, G5 and Gg as given 


below. 

G 7 = 2910 (3<0 (3.19a) 

G 7 = 2910 + 2700P - 343P 2 0 < p < 5.34 (3.19b) 

G 7 = 7560 |3 > 5.34 (3.19c) 

G 8 = -96900 J3<0 (3.19d) 

G g = - 988, 000 - 1.15 x 10 6 p + 89, 000|3 2 0 < p < 5.34 (3.19e) 

G 8 = -4.57 x 10 6 |3 > 5.34 . (3.190 

The value of (3 used to define G7 and Gg is lagged in time. 


The final model, called MMLPG, was developed using the equilibrium turbu- 
lent flows of Bradshaw and is defined by equations (3.14) through (3.19). The result- 
ing velocity-defect profiles for all three pressure gradient flows are shown in figure 17 
and exhibit good agreement with the experimental data, with the exception of the 
strong pressure gradient case. The calculations were performed on a CRAY Y-MP 
computer and the computational times are given in table 2. The strong pressure gradi- 








Table 2. Computational times for flat plate flows. 


(a) Zero Pressure Gradient 


Model 

Iterations 

sec./iter7 
grid point 

BLM 

MMLPG 

2000 

2000 

2.02x1 O' 5 
2.00x1 0* 5 


(b) Mild Pressure Gradient 


Model 

Iterations 

sec. /iter./ 
grid point 

BLM 

MMLPG 

3000 

3000 

2.14xl0" 5 

2.17zl0* 5 


(c) Strong Pressure Gradient 


Model 

Iterations 

sec ./iter./ 
grid point 

BLM 

MMLPG 

18,000 

10,000 

2.14xl0‘ 5 

1.96xl0* 5 
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ent case took considerably longer to reach convergence because the code had difficul- 
ties resolving oscillations induced at the upstream boundary, which was a fixed 
velocity profile as mentioned in section 3.3. 


3.5 Averaging for Multiple Boundaries 

If both walls in a given coordinate direction are solid surfaces, the turbulent 
mixing lengths are computed separately at each surface and then averaged. The 
Sajben diffuser, which is described in Chapter IV, has solid walls at the upper and 
lower vertical boundaries, and is a typical example of a geometry which would require 
averaging of the mixing length. The averaging formula of Appendix 3 (equation C.9), 
which was used to average the F wa i. c function in the Baldwin-Lomax model, is also 
used here to average the mixing length: 


Af I + ^2^2 
f,+f 2 


(3.20) 


If the lower and upper boundaries in the vertical direction, are solid surfaces, as in the 
Sajben diffuser, then l\ and 1 2 are the mixing lengths at the lower and upper bound- 
aries, respectively. The functions fj and f 2 are defined in equation C.10 of Appendix 
3. 



CHAPTER IV 


ADVERSE PRESSURE GRADIENT TEST CASES 

To evaluate MMLPG for some typical propulsion flows, a converging-diverg- 
ing duct, called the Sajben diffuser, was used. This duct is the diffuser portion of the 
inl et for a rocket/ramjet propulsion system; detailed experimental and computational 
data are available in the literature for flows both with and without external pressure 
pulse excitations. 30 ' 35 This study, however, dealt only with the unexcited flows. The 
geometry of the diffuser is given in figure 18: the throat height, H, is 44 mm; the 
entrance-to-throat ratio is 1.4, and the exit-to-throat ratio is 1.5. The grid, shown in 
figure 19, is the same as that used by references 13 and 34, and has 81 streamwise 
points and 51 vertical points. It was packed in the vertical direction near the walls in 
order to resolve the turbulent boundary layers and in the streamwise direction near the 
throat to resolve the normal shock. The reliability of this grid is discussed in Appen- 
dix 2. Three transonic flow cases were run. The flowfields were determined by 
setting R, the ratio of the exit static pressure to the inlet total pressure. The first case 
had a weak normal shock with R=0.82; the second case had subsonic flow throughout 
(no shock) with Rs=0.862, and the third case had a strong normal shock with R=0.72. 
The reference velocity, temperature, pressure and length used in Proteus were 4.72 m/ 
s, 292 K, 135 kPa, and .044 m respectively. These values match the values used in 
other numerical simulations of this flow. 13,3 ^’ 34 The initial conditions were zero 
velocity and constant temperature and pressure everywhere in the flowfield. Both 
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Figure 18. Illustration of the Sajben diffuser geometry. 



Figure 19. Computational grid for the Sajben diffuser. 
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cases were run using MMLPG and two implementations of the Baldwin-Lomax 
model. 

Two slightly different implementations of BLM were used because it was 
discovered, during the course of these calculations, that a slight change in the BLM 
coding, which occurred in the update of Proteus from Version 1.0 14 to 2.0, 13 can 
effect turbulent calculations. References 14 and 34 give results for the weak and no 
shock cases calculated with Version 1.0 of Proteus and using BLM. These calcula- 
tions were repeated in the current study using the current version of Proteus, Version 
2.0, and slightly different results were obtained. (These results will be presented later 
in this chapter.) Discrepancies in the results were caused by differences in the imple- 
mentation of the BLM model, more specifically, in the F wake function. Version 2.0 of 
Proteus computes F wake using equation (C.6) in Appendix 3, and version 1.0 of 
Proteus computes F wake as 


Fw.ke = min<| 


v F 


max* max 


2 y„ 


c wk v diff 


(4.1) 


which is the formulation stated in the original paper by Baldwin and Lomax. 4 The 
BLM implementations using equations (4.1) and (C.6) will be referred to as BLM1 
and BLM2, respectively. 

4.1 Weak Shock Case 

The weak shock case was used as an example case in the Proteus User’s 
Manual, 13,14 and therefore was run first in order to gain familiarity with running this 
type of flow. It was computed as described in reference 13: First the exit pressure 
was gradually reduced to R = 0.1338 to establish supersonic flow throughout the 
diffuser; then it was gradually raised to R = 0.82, the desired ratio to establish the 
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weak normal shock, and iterated until the solution was no longer changing apprecia- 
bly with time. A plot of the static pressure on the top wall at two locations, one 
upstream and one downstream of the normal shock, as the solution progresses is given 
in figure 20. This indicates that pressure reaches a steady state level, which, for 
practical engineering purposes, can be considered a converged solution. A closer 
examination of the results indicates that the solution oscillates slightly about a mean 
steady level. This may be caused by inherent unsteadiness in the flow; Salmon et al 30 
mention that very low-amplitude, self-sustaining oscillations were observed experi- 
mentally. It is more likely that the oscillations present in this calculation are numeri- 
cal in nature, which is common for flows with shock waves. The oscillations 
originating at the shock may not be damped out by the artificial viscosity and there- 
fore tend to migrate upstream. The artificial viscosity used in Proteus to calculate this 
flow was second- and fourth-order explicit, both using the nonlinear coefficient model 
of Jameson et al ; 19 the respective smoothing coefficients are K 2 and K 4 , as given in 
Appendix 2. For the entire calculation, K 2 was set to 0.1; K 4 was set to .005 for the 
first 6000 iterations, while the exit pressure was changing, and decreased to .0004 for 
the remaining 3000 iterations, which were at a constant exit pressure. More details 
about the effects of the artificial viscosity on this solution are included in Appendix 2. 

The static pressure distribution on the top and bottom walls is given in figure 
21. The small discrepancies just downstream of the shock are due to insufficient 
streamwise grid distribution. The shock location on the upper wall and the shock 
Mach number at the edge of the upper wall boundary layer are given in table 3. Both 
MMLPG and BLM2 do a good job of predicting the pressure distribution on the wall 
and the location of the shock. Each case was run for 9000 iterations and calculations 
made using MMLPG, BLM1 and BLM2 required 3.44 x 10 ' 5 sec/iteration/grid point, 
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(a) MMLPG 



(b) BLM2 


Figure 20. Static pressure history at two locations on the top wall: just upstream and 
just downstream of the normal shock. 
















Table 3. Shock location and Mach number, weak shock case. 


Turbulence 

Model 

Shock 

Mach 

Number 

Shock 

Location 

(x/H) 

MMLPG 

1.233 

1.57 

BLM1 

1.309 

1.73 

BLM2 

1.228 

1.49 

Experiment 30 

1.235 

1.41 
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3.90 x 10 -5 sec/iteration/grid point and 3.36 x 10“ 5 sec/iteration/grid point, respec- 
tively on the CRAY Y-MP computer. 

4.2 No Shock Case 

The second diffuser test case did not have a normal shock wave. To compute 
this case, the exit pressure was gradually lowered to R=0.862 then iterated until the 
solution stopped changing. A steady state solution was reached with subsonic flow 
throughout the entire diffuser. The Proteus default artificial viscosity, which uses the 
constant coefficient model of Steger 18 with both fourth-order explicit and second- 
order implicit artificial viscosity, was used; the smoothing coefficients, e ^ and £/ 
(defined in Appendix 2), had values of 1.0 and 2.0, respectively. 

The static pressure distribution on the top and bottom walls is given in figure 
22 and shows that MMLPG is clearly better than BLM1 and BLM2 at predicting the 
pressure distribution, though it still predicts a larger pressure drop than that given by 
the experimental data. The MMLPG results are similar to the calculations of Hseih et 
al 33 who attributed the lower throat pressure to the fact that the experiment was highly 
sensitive to small perturbations in exit pressure. The maximum Mach numbers in the 
diffuser are given in table 4. Though no experimental data is available to compare 
these values, the MMLPG results are in best agreement with the calculations of 
Georgiadis 35 who used the PARC Navier-Stokes code 36 for the same geometry. Each 
case was run for 9000 iterations and calculations using MMLPG, BLM1 and BLM2 
required 3.45x1 O' 5 sec/iteration/grid point, 2.90xl0" 5 sec/iteration/grid point and 
3.36xl0" 5 sec/iteration/grid point, respectively, on the CRAY Y-MP computer. 









Table 4. Maximum Mach number, no shock case. 



Maximum 

Turbulence 

Mach 

Model 

Number 

MMLPG 

0.881 

BLM1 

0.972 

BLM2 

0.976 




49 


4.3 Strong Shock Case 

The final diffuser flow computed was the case with a strong normal shock 
positioned in the throat. This case was run in a manner similar to the weak shock 
case: First the exit pressure was gradually lowered to R=0.1338 to achieve supersonic 
flow throughout the diffuser; then it was gradually raised to R=0.72 to establish the 
strong normal shock in the throat region, then iterated there until the solution stopped 
changing appreciably with time. Proteus was run in both steady and unsteady modes 
to try to simulate the experimentally observed self-excited oscillations of 217 Hz. 30 
Unsteady mode in Proteus is achieved by calculating a global time step whereas 
steady mode uses a local time step to speed up the computation. Neither calculation 
simulated the experimentally observed oscillatory behavior, but instead produced very 
small numerical oscillations in the flow properties. (The artificial viscosity used for 
the strong shock calculations was the same as that used for the weak shock calcula- 
tion.) Figure 23 shows the static pressure on the top wall at the experimental shock 
location and illustrates the behavior of these small oscillations; the calculation shown 
was run in unsteady mode with MMLPG. 

The static pressure on the top and bottom walls are presented in figure 24 and 
the shock location and Mach number at the edge of the top wall boundary layer are 
given in table 5. Both BLM2 and MMLPG predicted the shock location too far 
downstream, while BLM1 predicted it too far upstream. The experiment predicted a 
region of separation on the top wall just downstream of the shock with the flow 
reattaching at x/H=6.0. MMLPG predicted a very small region of separation on the 
top wall which reattached at x/H=3.6. BLM2 predicted very small regions of separa- 
tion on both the top and bottom walls which reattached at x/H=3.8 and x/H=6.2, 
respectively, and BLM1 predicted a separation along the bottom wall that did not 
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Figure 23. Shock static pressure on top wall for the Sajben diffuser, strong shock 
case. 
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Table 5. Shock location and Mach number, strong shock case 


Turbulence 

Model 

Shock 

Mach 

Number 

Shock 

Location 

(x/H) 

MMLPG 

1.626 

3.13 

BLM1 

1.411 

2.11 

BLM2 

1.665 

2.90 

Experiment 

1.353 

1.98 
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reattach. The separated behavior is illustrated in figure 25, which gives the velocity 
profiles at four locations downstream of the shock. These peculiar results can be 
attributed to the behavior of the turbulent viscosity as shown in figure 26. MMLPG 
and BLM2 compute much higher values of |X t than does BLM1; this is due to the large 
increase in vorticity downstream of the shock. In BLM1, as given in equation (4.1), 
the second formula results in F wake being inversely proportional to the vorticity 
magnitude, and this helps to reduce the value of p t for this model. Each case was run 
for 10,000 iterations and the steady calculations using MMLPG, BLM1 and BLM2 
required 3.54xl0' 5 sec/iteration/grid point, 3.83xlO' 5 sec/iteration/grid point, and 
3. 82x1 O' 5 sec/iteration/grid point, respectively, on the CRAY Y-MP computer. 




u(m/s) 


u(m/s) 


(c) x/H = 6.34 


(d) x/H = 7.49 


Figure 25. Velocity profiles for the strong shock case 
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Figure 26. Turbulent viscosity ratio, jLX t /jJ., for the Sajben diffuser, strong shock case. 



CHAPTER V 


SUMMARY AND CONCLUSIONS 

The objective of this work was to modify the MML algebraic turbulence 
model to increase its range of applicability to include zero and adverse pressure gradi- 
ent boundary layer flows. To accomplish this objective, modifications were made 
based on experimental data for zero and adverse pressure gradient flat plate flows. 

The resulting model, called MMLPG, successfully predicted the flat plate flows and 
also successfully predicted three transonic diffuser flows. This indicates that the 
objective of this work has been met. 

In order to provide meaningful solutions for turbulent flows, CFD codes 
require good turbulence models. Algebraic models are the simplest and the most 
computationally inexpensive of turbulence models, and so were chosen as the focus of 
this study. Proteus, which was used to make all of the calculations in this work, is a 
Reynolds averaged Navier-Stokes code for aerospace propulsion flows and contains 
the Baldwin-Lomax algebraic turbulence model as a default. The Baldwin-Lomax 
model is known to have problems calculating certain flowfields, namely flows with 
strong pressure gradients and large separated regions. A promising newer model, the 
MML model, produced significantly better results than the Baldwin-Lomax model for 
separated airfoil flows, 5 but it was not evaluated for other types of flows. The objec- 
tive of the current work was to continue development of MML to increase its range of 
applicability to include boundary layer flows with zero and adverse pressure gradi- 
ents. 
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To accomplish this objective, MML was installed in Proteus and first used to 
calculate zero pressure gradient flow over a flat plate. These results indicated that the 
original MML model was not allowing for the proper boundary layer growth with 
increasing Reynolds number. To remedy this behavior, a relationship was found for 
the saturation length scale as a function of the local skin friction coefficient. This 
modified version of MML allowed for the proper boundary layer growth and therefore 
produced the correct velocity-defect profiles. Next, MML was modified to calculate 
adverse pressure gradient flows using the experimental data of Bradshaw for zero, 
mild and strong adverse pressure gradient flows. These modifications were combined 
into one generalized model, called MMLPG. This new model accurately predicted the 
zero and adverse pressure gradient flows and exhibited better agreement with experi- 
mental data than the Baldwin-Lomax model. 

To more thoroughly evaluate MMLPG for other adverse pressure gradient 
flows, this model was also used to calculate three transonic diffuser flow cases: flow 
with a weak shock, flow with no shock, and flow with a strong shock. These are flows 
typically encountered in aerospace propulsion applications. The MMLPG results 
were compared with results calculated using two slightly different implementations of 
the Baldwin-Lomax model, referred to as BLM1 and BLM2. The differences in the 
two models arise from the calculation of the F wa k e function, as discussed in Chapter 
IV. 

For the weak shock case, MMLPG and BLM2 did equally well in predicting 
the shock Mach Number and location, and also in predicting the static pressure distri- 
bution on the top and bottom walls. However, BLM1 over-predicted the shock Mach 
number and location and did not match the wall static pressures in the throat region of 
the diffuser. For the no shock case, MMLPG was significantly better than either of the 
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Baldwin-Lomax models at predicting the static pressures on the walls and MMLPG 
also predicted a maximum Mach number that was in close agreement with the results 
of Georgiadis. 35 

The strong shock diffuser flow was beyond the capabilities of all three 
models. Both MMLPG and BLM2 over-predicted the shock Mach number and 
location, as well as the pressure distribution on the walls, while the BLM1 results 
were in reasonable agreement for these quantities. None of the models were able to 
correctly predict the shock-induced separation on the top wall; in fact, BLM1 
predicted no separation at all on the top wall yet predicted a very large, unattached 
separation on the bottom wall. The primary reason the MMLPG and BLM2 results 
differ greatly from the BLM1 results is explained by the turbulent viscosity values. 
MMLPG and BLM2 gave maximum turbulent viscosity values of 23,000 and 21,000 
times the molecular viscosity, respectively, while BLM1 gave a maximum turbulent 
viscosity of only 4,600 times the molecular viscosity. The turbulent viscosity is 
proportional to the vorticity, which becomes very large just downstream of the normal 
shock, however, the F wake function used in BLM1 is inversely proportional to the 
vorticity resulting in a lower overall turbulent viscosity. The poor performance of all 
of the models for this case can also be attributed to the fact that all of the models are 
equilibrium turbulence models being used to calculate a flow which is clearly 
nonequilibrium. The poor performance of MMLPG for the strong shock case is also 
explained by the derivation of the model, which is based on experimental data for P 
values between 0 and 5, while this flow encountered (5 values as high as 12,000. 
Considering that MMLPG did well at calculating the less severe no shock and weak 
shock flows indicates that it is most likely applicable to other propulsion flows at 
these less severe types of conditions. 
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The flat plate and transonic diffuser results indicate that the modified version 
of the MML model, MMLPG, is capable of accurately predicting turbulent flows with 
and without adverse pressure gradients. Thus, the objective of this work, which was 
to evaluate the original MML model and modify it to increase its range of applicabil- 
ity to include adverse pressure gradient flows, has been met. Future work should 
include continued validation of the model for these types of flows as well as continued 
development of the model to better account for stronger adverse pressure gradient 
flows both with and without separation. 
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APPENDIX 1 


GOVERNING EQUATIONS OF PROTEUS 

The governing equations of Proteus are the compressible Navier-Stokes 
equations. The equations given below are taken from reference 13, but may also be 
found in several references. 2,3,38 Since the two-dimensional/axisymmetric version of 
the Proteus code was used for all calculations discussed in this work, the two-dimen- 
sional, planar equations are given here. (For the axisymmetric version of the 
equations, which are somewhat more complex, consult reference 13.) 

1 . Cartesian Coordinates 

In Cartesian coordinates, the continuity, x-momentum, y-momentum and 
energy equations are written in strong conservation law form, using vector notation: 


where 


3Q 3E 3F _ 3F k 

dt 8x dy dx dy 



(A.1) 


(A.2) 
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E = 


pu 

pu 2 + p 
puv 

(E r + p) u 


(A.2a) 


pv 
puv 
pv 2 + p 
(E r +p) v 


E v 


0 


'xy 

UX + VX ■ 

_ v xx xy 


(A.2b) 


(A.2c) 


F„ = 


0 


xy 


yy 


^y+^yy-q^ 


(A.2d) 


In equation (A.2), the dependent variables are p, pu, pv and E r , the inviscid flux 
vectors are E and F, and the viscous flux vectors are Ey and Fy. The normal and 
shear stresses, and the heat .flux are given by 


_ 9u » ,9u , dv s 
x « “ 2 ^ +X,( 0^ 3y } 


_ 3v , . ,3u 3v. 
x yy " 2 ^ +X( 3x 9y } 


8u dv . 
x *y “ M ’ ( d^ + 9x ) 


(A.3) 
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,3T 

dx 


q y = -k 


3T 

dy 


2. Equation of State 

In addition to the above equations, an equation of state is needed to relate the 
pressure to the dependent variables. Proteus contains the equation of state for an ideal 
gas: 


p = pRT 


(A. 4) 


where R is the ideal gas constant. For a calorically perfect gas, this is equivalent to 

P = (Y-l) [E r -ip(u 2 + v 2 )] (A.5) 


3. Generalized Grid Transformation 

The governing equations in section 1.0 are written in Cartesian coordinates 
and are not well-suited for nonrectangular geometries and grids with unequal 
spacing. 13, 39 To overcome these difficulties, the following generalized grid transfor- 
mation is used to transform the governing equations from physical (x, y) coordinates 
to rectangular orthogonal computational (£, T|) coordinates. 


S = 5(x,y) 

11 = r| (x, y) 


(A.6) 


The resulting spatial computational domain is square, and has uniform spacing. The 
chain rule is used to transform the partial derivatives in the Cartesian form of the 
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governing equations, (A.1)-(A.3); for details, refer to Towne et al. 
equations are written as 

The transformed 

3Q BE BE _ 3E V 9F y 
8t + 8^ + dr[ - 9^ 9t| 

(A.7) 

where 



(A.8) 

E= I(E^ + F^ + Q^ t ) 

(A.8a) 

F = j(Eq x + Fq y + QTl t ) 

(A.8b) 

Ev = -j (Ev^ x + F v ^ y ) 

(A. 8c) 

F v = j(E v q x + F v n y ) 

(A.8d) 

and J is the Jacobian of the transformation. 


7= ^Tl y -^ y Tl x 

(A.9) 


4. Time Differencing 

The generalized time differencing scheme of Beam & Warming 16 is used to 
approximate the time derivative in equation (A.7): 
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3Q AQ - 8, 3(AQ") r 1 ASi- 

at At i+e 2 at i+e 2 at i+e 2 At 

+ o[(0 1 -i-0 2 )At, (At 2 )] 


where AQ n = q " +1 _ q" an( j the superscripts n and n+1 denote the known and 
unknown time levels, respectively. 

The parameters ©i and 02 determine the type of time differencing used. Table 
6 summarizes the available schemes. 


Table 6. Time differencing schemes in Proteus 


0i 

0 2 

Method 

Truncation Error 

1 

0 

Euler implicit 

' O(At) 

1/2 

0 

Trapezoidal implicit 

O(At) 2 

1 

1/2 

3 -point backward implicit 

O(At) 2 


The Euler implicit method is recommended for steady flows and the 3-point backward 
implicit method is recommended for unsteady flows. 

5. Space Differencing 

Spacial first derivatives in the £, direction are approximated using the follow- 
ing second-order central difference formula. 


s 5 A _ /i+l.j /i-l.j 

( 35V 2A^ 


(A. 11) 


The computational grid spacing, A^, is constant. A similar formula is used for first 
derivatives in the rj direction. 




APPENDIX 2 

ARTIFICIAL VISCOSITY AND GRID CONVERGENCE 


High frequency nonlinear instabilities can appear as the Proteus solution 
develops. For example, physical phenomena, such as shock waves, can cause instabil- 
ities when they are captured by the finite difference algorithm. In addition, high 
Reynolds number flows may have oscillations resulting from the odd-even decoupling 
inh erent in the use of central spatial differencing of the convective terms. Artificial 
viscosity may be used to suppress these oscillations. The two artificial viscosity 
models in Proteus are the constant coefficient model of Steger 18 and the nonlinear 
coefficient model of Jameson et al. 19 The implementation of these models in general- 
ized nonorthogonal coordinates was taken from Pulliam. 40 

1 . Constant Coefficient Model 

The constant coefficient model uses a combination of explicit and implicit 
smoothing. The standard explicit artificial viscosity uses fourth-order differences. 
Second-order explicit artificial viscosity, which provides more smoothing, is also 
available in Proteus, however it is rarely used because it introduces a large error. The 
implicit smoothing is second order and is used to extend the linear stability bound of 
the fourth-order explicit smoothing. 

The explicit artificial viscosity is implemented in the Proteus alternating 
direction implicit (ADI) algorithm 15 by adding the following terms to the right-hand 
side source term for the first ADI sweep. 
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ef At 


(V^Q + V^Q) 


e< 4) At 

~~T~ 


[(V^) a Q+-(VA) a Q] 


(B.l) 


where t E ^ and e^ 4) are the second- and fourth-order explicit artificial viscosity 
coefficients. The symbols V and A are the backward and forward first difference 
operators for the first ADI sweep such that 


V^Qi = Qi-Qi_i 
A § Qi = Q i+ i-Qi 
V^Q i = Q i+1 -2Q i + Q i _ 1 
( V^) 2 Qs = Qi+2 _4 Qi+i + 6Qi - 4Q ; . , + Qi-2 


(B.2) 


Similar formulas are used in the tj direction. 

The implicit artificial viscosity is implemented by adding the following terms 
to the left-hand side of the governing equation. 

e,At . 

— — [V 5 A 5 (/AQ )] (B.3a) 


£,At 

~T 


[V„A„(7AQ n )] 


(B.3b) 


Equation (B.3a) is added for the first ADI sweep and equation (B. 3b) is added for the 
second ADI sweep. The constant £/ is the implicit artificial viscosity coefficient. 

The optimum values of the coefficients Ej/ 4 -* and £/ vary from problem 
to problem. They should be small so as not to corrupt the physical solution, yet large 
enough to damp any instabilities. The Proteus User’s Guide ^ recommends starting 
values of E^ 4 -* =1.0, ££^=1.0 and £/= 2.0. 
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2. Nonlinear Coefficient Model 

The nonlinear coefficient artificial viscosity is explicit and contains second 
and fourth-order differences. The following terms are added to the right-hand side of 
the governing equations. 


v,{ 


(?) + ( 7> J (e “ A « Q " 

+ { [t,., + 7,] 


The difference operator A^V^A-Q; is defined by 

A|V|A^Qi = Qi + 2 _ 3 Qi+i + Q1-1 


and the expression \|/ is given by 


¥ = ¥* + ¥ y 


where \|/ x and \|/ are spectral radii defined by 


(B.4) 


(B.5) 


(B.6) 


¥ v : 


IUI+a M+% 

M 

\V\+ahl + T\ 2 y 
Vt| 


(B.7) 


The second- and fourth- order nonlinear artificial viscosity coefficients are a 
function of the pressure field. In the £ direction, they are given by 


(ef). = K 2 At max(c i+1 , Oj, a^,) 


(ej 4) ) . = max [0, K 4 At - (e| 2) ) .] 


(B.8a) 

(B.8b) 


where 
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<Jj = 


p i+ i + 2 pi + pi-i 


(B.9) 


Similar formulas are used in the tj direction. 

The parameter a is a pressure gradient scaling parameter which increases the 
amount of second-order smoothing relative to fourth-order near shock waves. The 
parameters k 2 and K 4 are user-specified constants. As with the constant coefficient 
model, the optimum values of k 2 and K 4 are problem dependent. Typical values range 
from k 4 =0.005 and k 2 =0.01 for flows with no shocks, to k 4 =0.0004 and k 2 =0.1 for 
flows with shocks . 13 Pulliam gives k 2 =0.25 and k 4 =0.01 as typical values for an 
Euler analysis . 13, 40 

3. Comments on Artificial Viscosity 

As previously mentioned, artificial viscosity is generally used to minimize 
oscillations which occur when computing high Reynolds number flows and flows with 
shock waves. Since the artificial viscosity terms do not represent anything physical, 
the coefficients should be as small as possible so as not to corrupt the solution, yet 
large enough to damp the nonphysical instabilities. Optimum values of the artificial 
viscosity coefficients vary from problem to problem; the coefficients used to calculate 
the flows presented in Chapters III and IV were selected based on values used for 
similar cases, as given in the Proteus User’s Manual . 13 Some representative test cases 
were evaluated to insure that the chosen artificial viscosity did not corrupt the physi- 
cal characteristics of the flow. 

The flat plate flows presented in Chapter III were run using the constant 
coefficient model with e £ ( 4 ) =1.0, e £ ( 2 ) =0.0 and e 7 =2.0. For these flows, it was possi- 
ble to run Proteus with zero artificial viscosity, however the solutions took two to four 
times longer to “converge,” or reach a point where the solution stopped changing 
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appreciably with time. Upon close examination, these solutions did not agree as 
closely with experimental data as the solutions computed using artificial viscosity. 

For the diffuser flows computed in Chapter IV, the artificial viscosity effects 
were examined for the weak shock case. As mentioned in Chapter IV, the nonlinear 
coefficient model was used. A value of k 2 =0.1 was used for the entire calculation, 
with K4=0.005 while the exit pressure was changing (i.e, for the first 6000 iterations), 
and k 4 =0.0004 for the remaining 3000 iterations, which were at a constant exit 
pressure. It was not possible to compute this flow without artificial viscosity, so the 
effects of doubling and halving the smoothing coefficients was examined. The static 
pressure distribution on the top and bottom walls for this comparison (computed using 
MMLPG) is given in figure 27. The solution computed using half of the original 
artificial viscosity was nearly identical to the original solution, indicating that the 
originally chosen artificial viscosity is reasonable for this flow. Doubling the artificial 
viscosity gave a less desirable result in that the normal shock was smeared over a 
greater number of grid points. 

4. Grid Convergence 

Grid convergence is an important factor in the accuracy of a CFD calculation. 
The grids used to make the flat plate and transonic diffuser calculations were assessed 
to insure their grid independence. For the zero pressure gradient flat plate calcula- 
tions, a 101x101 grid was initially chosen. The size of this grid was systematically 
reduced in each direction in order to find the coarsest grid that would give a solution 
which would not change if additional grid points were added. The 51x51 grid 
described in Chapter III was chosen based on this study. 

The grid used to make the transonic diffuser calculations had been used previ- 
ously by others, 13 ’ 34 so it is probable that this grid gives a reliable solution. As an 



'5 



Figure 27. Comparison of weak shock static pressure distributions, computed using 
MMLPG and three different amounts of artificial viscosity. 
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added check, the number of grid points in each direction was doubled, and the result- 
ing 162x101 grid was used to compute the no shock flow using MMLPG. A compari- 
son of these results with the results obtained using the 81x51 grid of Chapter IV is 
given in figure 28 and indicates that the 81x51 grid is reliable. 
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(a) Top Wail 



(b) Bottom Wall 


Figure 28. Comparison of no shock static pressure distributions, computed using 
MMLPG and two different grids. 





APPENDIX 3 

THE BALDWIN-LOMAX TURBULENCE MODEL 


A generalized version of the Baldwin-Lomax algebraic turbulence model is 
available in Proteus. 13 As mentioned in chapter in Chapter II, the turbulent shear and 
normal stresses and the turbulent heat flux are modeled using the Boussinesq 
approach, where the effective viscosity is defined as JX eff = p + p t , the second coeffi- 
cient of viscosity is defined as X eff — X + A, t , and the effective thermal conductivity 

coefficient is defined as k eff = k + k t . 

For wall bounded flows, the Baldwin-Lomax model is a two-layer model: 


= 


dinner’ ? * 


outer ’ y ^crossover 


(C.l) 


where y croS sover is smallest value of y at which the inner and outer region values of p t 
are equal. For free turbulent flows, jx t = ({i t ) outer - 


1. Inner Region 

The inner region turbulent viscosity is computed from 


(ft.). = p/ 2 !®! 

inner 1 1 1 


(C.2) 


where l is the mixing length given by 
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( 


l = Ky 


l 1 



(C.3) 


The quantity |<n| is the magnitude of the total vorticity, defined for two-dimensional 
planar flow as 


| ©| 


dv 9u 
9x 9y 


(C.4) 


2. Outer Region 

In the outer region, the turbulent viscosity is given by 

(M- t ) outer = ^C cp pF Kleb F wake (C.5) 

where K is the Clauser constant, taken as 0.0168 and C C p is a constant taken as 1.6. 
The quantity F wake is computed from 


p — 

■•■wake 


y maxF n 


n v 2 y max 

•^wk v diffp 

A may 


for wall bounded flows 


for free turbulent flows 


(C.6) 


where the constant C w j. is 0.25 and 


v d ,„ = lv|»„-|vL,„ 

where V is the total velocity vector. The quantity F max is the maximum value of 


F(y) = 


y|©| 


1-e A+ 
V 


iox wall bounded flows 


(C.7) 


L y|co| . 


for free turbulent flows 
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and y max is the value of y corresponding to F max . FjQ e b is the Klebanoff intermittency 
factor which accounts for the experimentally observed phenomenon that as the free 
stream is approached, the fraction of time the flow is turbulent decreases. It is given 
by 


Fjcieb “ 


1-5.5 (5^) 


i-i 


where CjQ e b is a constant taken as 0.3. 


(C.8) 


3. Multiple Boundaries 

If both walls in a given coordinate direction are solid surfaces, the turbulence 
model equations are applied separately at each surface and then averaged. The two 
outer regions overlap, and it assumed that the two inner regions do not overlap. The 
averaging procedure deals with the F^^g function. For example, in the vertical direc- 
tion, if the upper and lower boundaries are both solid surfaces, the two values of F^^g 
at a particular streamwise station are combined using the following averaging 
formula: 


wake 


(F W ake) i f l + (Fwake) 2 f 2 

fl + f 2 


(C.9) 


The quantities (F wake ) , and (F wake ) 2 are the separate values computed at the lower 
and upper boundaries using equation (C.6). The functions fj and f 2 are defined by 


fi« ( 


2D," 

yi } 


f 2 = ( 


2D," 

y 2 


(C.10) 
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The constant n is set equal to 2.0, yj and y 2 are the normal distances to the bottom and 
top surfaces, respectively, and D t and D 2 are the normal distances from the two 
surfaces to the location of | vL ax . In addition, the y/y max value used in equation (C.8) 
for FjQgb is computed for both surfaces and the minimum value is used. These values 
of F K i eb and F wake are then used in equation (C.5) to compute (H t ) outer . 

4. Turbulent Values of X and k 

The turbulent second coefficient of viscosity is defined as 



The turbulent thermal conductivity coefficient is defined using the Reynolds 
analogy as 

k, = ^ (C.12) 

Pr t 

and c p is the specific heat at constant pressure and Pr t is the turbulent Prandtl number. 
In Proteus, the turbulent Prandtl number may be equal to a constant or computed using 
the empirical formula of Wassel and Catton. 37 For the cases described herin, Pr t was 
constant with the Proteus default value of 0.91. 
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